arXiv:hep-ph/0503038v4 18 Aug 2005 


Can Thermalization in Heavy Ion Collisions be 
Described by QCD Diagrams? 

Yuri V. Kovchegov* 


Department of Physics, The Ohio State University 
Columbus, OH f 3210 

March, 2005 


Abstract 

The onset of thermalization in heavy ion collisions in the weak coupling framework 
can be viewed as a transition from the initial state Color Glass Condensate dynamics, 
characterized by the energy density scaling like e ~ 1/r with r the proper time, to the 
hydrodynamics-driven expansion of the quark-gluon plasma with e ~ (or higher 

power of 1/r for the boost non-invariant case). We argue that, at any order of the 
perturbative expansion in the QCD coupling constant, the gluon field generated in an 
ultrarelativistic heavy ion collision leads to energy density scaling as e ~ 1/r for late 
times r ^ l/Qs- Therefore it is likely that thermalization and hydrodynamic description 
of the gluon system produced in heavy ion collisions can not result from perturbative QCD 
diagrams at these late times. At earlier times with r ~ 1/Qs the subleading corrections 
to e in 1/r expansion (terms scaling like ~ l/r^"*"^ with A > 0) may become important 
possibly leading to hydrodynamic-like behavior of the gluon system. Still, we show that 
such corrections do not contribute to the particle production cross section, and are likely 
to be irrelevant for physical observables. We generalize our results by including massless 
quarks into the system. Thus, it appears that the apparent thermalization of quarks and 
gluons, leading to success of Bjorken hydrodynamics in describing heavy ion collisions at 
RHIC, can only be attributed to the non-perturbative QCD effects. 


e-mail: yuri@mps.ohio-state.edu 



1 Introduction 


Understanding how the system of quarks and gluons produced in ultrarelativistic heavy ion 
collisions evolves towards thermal equilibrium is one of the central questions in our theoretical 
understanding of nuclear collisions. On the one hand, there exists a strong experimental evi¬ 
dence for equilibration of the quark-gluon system produced in a heavy ion collision at RHIC. 
The evidence is based on the success of hydrodynamic models of the collisions P El 0 in E], 
indicating a collective behavior of the quark-gluon system, and on the discovery of jet quench¬ 
ing jni El IHl El EH EH 1121 uni eh eh lini im EHI, which demonstrated the presence of strong 
hnal state interactions. However, hydrodynamic models of the evolution of the quark gluon 
system only work well, especially for the elliptic flow V 2 PI. if equilibration occurs in fact at a 
very early time pin, t<0.5fm/c, a time whose smallness is difficult to reconcile with current 
dynamical pictures of equilibration [201 EH 1221 EHI EH EH EE| • 

On the other hand, a complete theoretical understanding of thermalization is still lacking. 
The success of saturation/Color Glass approach j2H EH EH EH EH E21 EH EH EH EH in de¬ 
scribing particle multiplicities PI in heavy ion collisions and particle spectra in deuteron-gold 

collisions |3HEHI2HliHI221llHlinilHI2HllHllHllHEHI^ (see also j02lEH) appears to 

indicate that saturation/Color Glass formalism is valid for the initial stages of heavy ion colli¬ 
sions at RHIC. Our understanding of the very early pre-equilibration stages of the collision and 
their description in terms of classical gluon helds has signihcantly advanced in the recent years 
ISnEHEHEHEHEHEHinH- Based on this saturation initial conditions, Baier, Mueller, Schiff 
and Son proposed the so-called “bottom-up” thermalization scenario EH in which multiple 
2 —> 2, 2 —> 3 and 3 —2 rescattering processes, the importance of which was originally under¬ 
lined in [02], would drive the system to thermal equilibration over the time scales of the order 
of To ~ While the estimates in PI were mostly parametric, the numerical value 

of this thermalization time appears to be much larger than 0.5 fm needed by hydrodynamic 
simulations PP. 

More recently it was argued by Arnold, Lenaghan and Moore that the “bottom-up” ther- 
malization scenario could be susceptible to plasma instabilities [HH En, which were advocated 
previously in [HH EH EH EH EH • Such instabilities might help facilitate the equilibration process 
making the thermalization time shorter than predicted by the “bottom-up” scenario. However, 
in En Arnold and Lenaghan proved a lower bound on the thermalization time, which happened 
to be surprisingly close to the “bottom-up” estimate. In PI it has been suggested that, while 
complete thermalization may not happen until later times, an isotropization of the produced 
particle distribution in momentum space may happen much faster, leading to generation of 
longitudinal pressure needed for hydrodynamic description to work. 

Here we take a different approach to the problem of thermalization. Thermalization could be 
thought of as a transition between the initial conditions, which are characterized by the energy 
density scaling like e ~ 1/r, and the Bjorken hydrodynamics, which, in case of the ideal gas 
equation of state has e ~ p. (Of course at realistic temperatures achieved in heavy ion 

collisions the power of 4/3 may become somewhat smaller: however, it is always greater than 
1 for hydrodynamic expansion.) Therefore it appears that corrections to the saturation/Golor 
Glass initial conditions [Sn EH EH EH EH EH would contribute towards modifying the e ~ 1/r 
scaling to some higher power. Thus one should be interested in Feynman diagrams which would 
bring in r-dependent corrections to e ~ 1/r scaling of the (classical) gluon helds in the initial 
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stages of the collisions. Unfortunately, after examining a number of diagrams, we noticed that 
while many of them introduce r-dependent corrections to the initial conditions, such corrections 
are subleading and small at large r and do not modify e ~ 1/r scaling at late times. After 
reaching this conclusion we have constructed a general argument proving that e ~ 1/r scaling 
always dominates at late times, both for classical helds and quantum corrections, which we are 
presenting here. 

The paper is structured in the following way. We begin in Section 2 by calculating the 
energy density of a lowest-order non-trivial classical gluon held from |HSj. As expected the 
energy density of the classical held scales as e ~ 1/r. We then continue in Section 3 by 
considering the most general case of boost-invariant gluon production, which is, indeed, not 
limited to classical helds. We argue that e ~ 1/r scaling persists to all orders in the coupling 
constant as shown in Eq. The argument is based on a simple observation (see Eq. (EH)) 
that r-dependent corrections to the classical gluon held may only come in through powers of 
gluon virtuality k‘^ in momentum space with each power of giving rise to a power of 1/r. In 
order for the on-mass shell amplitude (at = 0) to be non-singular only positive powers of /c^ 
are allowed: hence, the corrections come in only as inverse extra powers of r and are negligible 
at late times. In Section 3 we generalize our results to rapidity-dependent distributions. The 
e ~ 1/r scaling does not get modihed by rapidity-dependent corrections either (see Eq. (IH!?|l h 
Rapidity-dependent corrections come in as, for example, powers of fc+, which is one of light cone 
components of the gluon’s momentum. However, as could be seen from, say, Eq. (IFfI) , powers 
of do not modify the r-dependence of energy density. In Section 4 we argue that e 1/r 
scaling persists even when massless quarks are included in the problem. Therefore it appears 
that perturbative thermalization can not happen in heavy ion collisions. We conclude in Section 
5 by arguing that if perturbative thermalization is impossible, than the non-perturbative QCD 
effects must be responsible for the formation of quark-gluon plasma (QGP) at RHIC [Z^. We 
list the non-perturbative effects which we believe may be responsible for thermalization. 


2 Energy-Momentum Tensor of Classical Gluon Field 


We start by calculating the energy-momentum tensor of the lowest order gluon held produced 
in an ultrarelativistic heavy ion collision. This held has been found analytically in PI EH and 
the corresponding Feynman diagrams are depicted here in Fig. ^ The cross in Fig. ^denotes 
the space-time point in which we measure the gluon held. The gluon held in = 0 covariant 
gauge given by diagrams in Fig. Q can be written as p] 



(27r)"^ kd + ieko 




( 1 ) 



Figure 1: Lowest order (~ g^) gluon held produced in nuclear collisions. 
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( 2 ) 


where Cf^{k, q) is the Lipatov vertex [73j 

^'‘*'=•2) = + (3) 

The field of Eq. is given for a collision of a quark i in one of the nuclei having transverse 
coordinate Xj and light cone coordinate Xi- with a quark j in the other nucleus having transverse 
coordinate y. and the light cone coordinate yj^. The matrices {T^) and (7J) act in the color 
spaces of the quarks i and j correspondingly. Indeed we sum over all quark pairs in Eq. ©, 
with valence quarks i being in any of the Ai nucleons in the first nucleus and valence quarks j 
being in any one of the A 2 nucleons in the second nucleus. 

The energy-momentum tensor of a gluon field is given by 
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— _pa-y.p pav _|_ gpy ^pa Y 


We need to calculate averaged in the wave functions of both nuclei P?n IKK] 

1 




(4) 

(5) 


where the averaging implies integrating over all possible positions of quarks in the nucleons 
and nucleons in the nuclei, and taking traces (divided by NY) in the color spaces of the quarks 
[KTl IKK] . The averaging can be represented as 


(...) 


d'^Xi dYyj dxi- dyj+ 
Si± S 2 ± a- a+ 


]^Tr4Tr,[...] 


( 6 ) 


where S'l^ and S 2 ± are the cross sectional areas of the two nuclei, which we for simplicity 
assume to be cylindrical with the cylinder axis pointing in the beam [z-) direction. a_ and a+ 
are Lorentz-contracted nucleon sizes in the — and -|- directions correspondingly, which are very 
small, making averaging over Xi- and yj+ equivalent to just putting Xi- = 0 and yj+ = 0 jKK] . 

We are interested in calculating for the gluon field produced in a central nuclear collisions 
in the forward light cone. Since the gluon field of Eq. 0 is o{g^), we may only use it to compute 
o{g^) contribution to the energy-momentum tensor, for which we will need only the Abelian 
part of the field strength tensor To compute higher orders in one would also need 
higher orders in A“. Using the field (HD in the Abelian part of Eq. (0, performing the averaging 
defined in Eq. 0 and remembering that the multiplicity distribution given by the diagrams of 
Fig.E for gluons with transverse momentum k, rapidity y, located at impact parameter b, is 

dN A 1 A 2 1 kp 

d?kdyd% vr SipS 2 x^ A 
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we obtain after lengthy algebra (and dropping the averaging sign around 
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( 8 ) 


Here we used x± = (t±z)/\/2, r = — z‘^ = y/2x+x-, and kx = \k\. We also took advantage 

of the fact that the multiplicity distribution m is rapidity independent and, since should 
depend only on space-time coordinates, replaced momentum space rapidity y with the space- 
time rapidity rj = (1/2) ln(a:+/a:_). (At this point such substitution makes no difference: in 
Section lO we will show how this substitution is formally justihed in the rapidity-dependent 
case.) In arriving at Eq. (jHl) we have used the integral dehned in Eq. ea and the one given 
by Eq. dBH) in Appendix B with A = 0 along with 


d^q dvr ^ kr 

q\k-qy ^ 


(9) 


where A is some infrared cutoff. Eq. (jH)) is derived in the leading logarithmic approximation in 
In kx/A.. 

Eq. (jH)) gives us the energy-momentum tensor in the forward light cone of the lowest order 
gluon held from Fig. C] produced in a central collision of two identical nuclei. While it is written 
in a non-specihc form with regards to the order of the coupling constant g, we have proven Eq. 
(|H|) only at the order o{g^). 


3 Energy Density in the Boost-Invariant Approximation 


3.1 Region of Applicability 


Let us hrst consider the case of high energy heavy ion collisions, where the total rapidity interval 
is large enough to allow for eikonal approximation, but not large enough for quantum BFKL- 
type corrections ca to become important. This is the quasi-classical regime of McLerran- 
Venugopalan model isniiHiiiHa. To achieve it one needs the Bjorken x variable to be small 
enough such that na 


X < 


1 

2m]\fR 




( 10 ) 


which is the condition ensuring that coherent eikonal interactions are possible in the nuclear 
wave functions. For corresponding rapidities, Y = In l/x, the condition of Eq. (fUH) means that 


Y > InA^/^ 


( 11 ) 
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Remembering that McLerran-Venugopalan model corresponds to resummation of multiple 
rescatterings parameter ~ 1 we rewrite Eq. m as 

Y > In — ~ In—. (12) 

aj tts 

On the other hand, the boost invariant approximation is broken down by quantum evolution 
corrections, which, in the dominant leading logarithmic approximation, bring in powers of agY 
[I31ISS1IS11ISE|. Indeed these corrections are negligible when asY <1, such that 

r < —. (13) 

(y.g 

Eqs. m and m dehne the rapidity interval for nuclear collisions in which the boost invariant 
approximation, which we will consider in this Section, is valid. 


3.2 Most General Boost Invariant form of 

Similar to the Bjorken approach jl] we will consider a central collision of two very large nuclei, 
such that the gluon production is translationally invariant in the transverse direction. Dehning 
two four-vectors in terms of light-cone coordinates 

(14) 

and 

(15) 

we can write the most general energy-momentum tensor for the system as 


= A{t) + B{t) {u^Vy + UyVf,) + C{t) v^Vy + D{t) g^y, (16) 


where in our convention g^^y = diag{l, —1, —1, —1}. Here the parameters A, B, C, D in Eq. |Tn|) 
are functions of r only, since the large transverse extent and azimuthal cylindrical symmetry 
for central collisions of the nuclei allow us to neglect the transverse coordinate dependence, 
and the assumption of boost invariance makes functions A, B,C, D independent of space-time 
rapidity g = (1/2) ln(a;+/a;_). 

From Eq. (uni) we see that 


r++ = WT) + B(r) + C(r)](/t)' 


(17) 


and 

T.. = [A{T)-B(r) + C{r)](^'j\ (18) 

Due to -f — symmetry of the collision of two identical nuclei it should be possible to obtain 

T_from r++ after changing all -|- indices to — indices of all the relevant four-vectors in it. 

This condition, when applied to Eqs. ira and ca, demands that B{t) = 0. 
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Rewriting 


we obtain for 


the remaining non-zero functions A, C, D as 

A{t) = e(r) p(r), C(r) = p^ir) - p(r), and D{t) = 

the non-zero components of the energy momentum tensor 

T++ = [e(r)+P3(r)] , 

T— = [e(r)+p 3 (r)] , 

T+- = [e(r) -psM] = [e(r) -psM] 

Tij = ■5yp(r), 


-p(r) 


(19) 


( 20 ) 


where the indices i,j = 1,2 denote the transverse components of the tensor. Eq. (EHD gives us 
the most general boost-invariant energy-momentum tensor for a collision of two very large nuclei 
with the total rapidity interval satisfying conditions m and m allowing for a boost-invariant 
description of the gluon production. 

At 2 ) = 0 in the center-of-mass frame the energy-momentum tensor from Eq. (jUl) can be 
written as 




f e(^) 

0 

0 

0 

\ 

0 

p(r) 

0 

0 

0 

0 

p(r) 

0 


1 0 

0 

0 

P3{t] 

/ 


( 21 ) 


Now we can see the physical meaning of the parameterization introduced in Eq. (P|) : e is the 
energy density and pa is the longitudinal pressure along the beam axis (^(-direction), which in 
principle does not have to be equal to the transverse pressure p. Indeed, for the case of boost- 
invariant Bjorken hydrodynamics P, the two pressures are identical, Ps{t) = p{t). However, 
as one can show, for classical gluon helds generated in heavy ion collisions (HU |S3 EHl EH EH] , 
the longitudinal pressure component is zero at sufficiently late times, Psir) = 0, while e(r) = 

2p(r)^0 |Sn|. 

Applying the conservation of energy-momentum tensor condition 



= 0 

(22) 

to the tensor in Eq. (On|l we obtain 

de _ e + p3 
dr T 

(23) 

similar to Bjorken hydrodynamics p]. 




Eq. shows that if energy density scales with proper time as e ~ 1/r then the longitudinal 
pressure is zero, ps = 0. This is indeed the case for classical gluon production in the initial 
stages of the heavy ion collisions considered in Section 2 (see also [521) • To see this let us use 
the energy momentum tensor from Eq. (jHj) in Eq. (EnD to write 

' = i I {[^i(fcTr)l^ + |Jo(4Tr)l=} 
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P = 


TT 


d^k 


dN 


d^kdr] d?b 


kr [MkTT)Y 


Using the large-argument asymptotics of the Bessel functions we write 


drk 


dN 


T>l/(fcT> 


T 


d?k dr] d?b 


kT = 


1 dEj' 


T drjd?E 


(24) 

(25) 


which precisely agrees with the Bjorken energy density estimate PJ. Here we assumed that the 
gluon spectrum is characterized by some typical transverse momentum (/cr), such that large 
time asymptotics is dehned by r 3> l/(/cr)- (Strictly speaking such “typical” momentum for 
lowest order gluon held of Eq. (HD is the infrared cutoff A, but it would become the saturation 
scale Qs A once multiple rescatterings are included mm-) 

Similarly, using the large-argument asymptotics of the Bessel functions one can show that 


P3 


r>l/(A:T> 


0 


(26) 


in agreement with Eq. and Eq. (123D- Thus the large time asymptotics of the energy- 
momentum tensor of the lowest order classical gluon held is given by = diag{e,p, p, 0} with 
e = 2p and e given by Eq. (ESD Similar results were obtained in numerical simulations of the 
full classical gluon held including all orders in multiple rescatterings |59j . 

The onset of thermalization or isotropization of the system EDI should come with generation 
of the non-zero longitudinal pressure ps comparable to the transverse pressure p. In order for 
that to happen Eq. (121 necessarily requires the energy density to start scaling with r as 
e ~ where A is some positive number. In the case of ideal Bjorken hydrodynamics 

A = 1/3. 

Thus the process of thermalization in heavy ion collisions can be viewed as a transition 
from the e ~ 1/r scaling, characteristic of free-streaming classical helds fl25jl . to e rsj 
scaling. Below we are going to study whether such transition can result from Feynman diagram 
resummation. 


3.3 Can Boost-Invariant Bjorken Hydro Result from Feynman Di¬ 
agrams? 

Let us explore what kinds of energy-momentum tensor may result from Feynman diagram 
resummation. We will concentrate only on gluon helds, and later will generalize our conclusions 
to include quark helds as well. We will assume that the initial gluon held is given by the classical 
held of McLerran-Venugopalan model [2Ul I2n EH EH EH EZl 1^1^ - though our results would 

^It is interesting to point out that in approaching the asymptotics of Eqs. and (ESI) both e andp 3 oscillate, 
such that e/3 becomes temporarily comparable to ps at proper time r ~ 1/Qs- While the mathematical origin 
of these oscillations is clearly due to the Bessel functions in Eq. EH, their physical interpretation (if it exists) 
is presently unclear. 
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not depend much on this assumption^. In the saturation scenario, the gluon helds at the early 
times with r ~ 1/Qs are strong 


A 


a 


rsj 


0^ 

9 


(27) 


In calculating corresponding held strength tensor one would require both the Abelian and 
the non-Abelian parts of it. However, as classical helds and their energy density (123) decrease 
with proper time, for t>1/Qs the Abelian part of would dominate. This is also true 
for quantum corrections to classical helds. Therefore, in the following discussion we will hrst 
restrict ourselves to calculating the Abelian part of only, and will later show that inclusion 
of non-Abelian parts of would not change our argument. 

The most general gluon held generated through any-order Feynman diagrams in = 0 
covariant gauge can be written as 



e 

(27r)"^ -|- ieko 


■m), 


(28) 


where we are using the retarded regularization for the outgoing gluon propagator —i/k'^ to 
ensure causality. The function denotes the rest of the diagram (the truncated part). In 

general J“(A;) is an abbreviated notation for J“(fc; qi, Ai; q 2 , A 2 ;...), which depends on momenta 
qi of extra gluons (or quarks) in the hnal state and on their polarizations (helicities) A*. In 
case of the classical gluon held there are no extra particles in the hnal state and momenta gds 
do not enter the expression ()28fl . (In fact, the possible particles in the hnal state for classical 
helds can be removed by using retarded regularization of gluon propagators 1ZS|.) Quantum 
corrections to the classical gluon held would inevitably bring in extra hnal state particles. 
The resulting “quantum” held A“(a;) in Eq. (OH|l would also depend on momenta gp again we 
suppress this dependence in the notation. Indeed, once quantum corrections are included there 
is no dominant gluon held anymore: in that sense the gluon held in Eq. (|2HD is not really a held, 
but more like a scattering amplitude (located on one side of the cut), with one {k) of the many 
outgoing particle lines (gds) being oh-mass shell with its propagator ending at a space-time 
point x^. 

Substituting the held from Eq. into the expression for the energy-momentum tensor 

(13 

{Tn = + 


and keeping only the Abelian parts of E“^’s we obtain 


T = 


d^kd^k' 


-.—ik-x—ik'-x 


(27r)® (A;2 -p ieko) + Fk'o] 


Ik^r'-ik) - k'-Mk)] [k'M') - k'j;(k')] 


+ j ik'r-ik) - k‘’r‘-(k}] ik;j;(k') - k',j‘^(k'}]), (29) 

^ There is a common misconception in the community that in McLerran-Venugopalan model one assumes 
that y = T]: while this assumption was made in the original works on the subject [HiEni, it is actually not 
necessary, with all the results of McLerran-Venugopalan model easily derivable without making any assumptions 
on correlations between 77 and y (see ED). 







where the brackets (...) are dehned by Eq. and now also inclnde integration over all momenta 
gds. The glnon held in Eq. (OH|) is generated by the color sonrces in two colliding nnclei, which 
are modeled by valence qnarks, jnst like in McLerran-Venngopalan model 1301 . Of conrse, the 
resnlting glnon held from Eq. (I2HD is not necessarily classical, it inclndes extra qnark and glnon 
emissions as well as loops, jnst like any prodnction diagram with incoming valence qnarks of 
the nnclei providing the initial condition for the scattering process. 

Since performing the transverse averaging over a very large nnclens in the brackets on the 
right hand side of Eq. (Ei pnts k = — fc', which resnlts from transverse translational invariance 
of glnon prodnction, we can rewrite it as 


T 


d^k d^k' ^-ik-x-ik'-x 

(27r)® + ieko) {k'"^ + iek'^) 


[k,r'>(k) - Fj“(fc)] [kU';(k') - k;j;{k')] 


+ Jft.. - k-J-'^ik)] lkX(k') - k'j;(k')]^y^S{k + kf), (30) 

where the donble brackets {{■■■)) denote now the color averaging, integration over g^’s, snm- 
mation over nncleons in the nnclei and averaging over longitndinal and remaining transverse 
coordinates. S± is the cross sectional area of the nnclei which we assnme to be identical. 

Let ns dehne the following correlation fnnction 


= ((J„(fc) Mk'))) 


(31) 


Using the covariant gange condition = 0, which translates into k^ J'^{k) = k'^ J^{k') = 0, 
along with the <->■ (and k'j^ <->■ k'_) symmetry of the collision, we can derive the following 
relations between different components of 


D+i — 

! Dii 

2 k- ^ 

D.^ = 

/ Dii 
2k+ 


K 


K 

Di+ = 

^ D 

2fc'_ 

A- = 

^ D 

2 A ' 


D++ 

D— 



k+ k'+ k_ k'_ ’ 


(32) 


where the Latin indices i,j = 1,2 indicate the transverse components of the correlators and 
two lowercase repeated Latin indices indicate contraction over that index. 

We are interested in calcnlating the energy density e given by (see Eq. ®) 


Using Eq. ® we write 


and 



(33) 

(34) 

(35) 
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Using Eq. ()28|1 together with relations from Eq. (jnS) in Eqs. (jHlI) and (jSSl) we obtain 


T++ = 


d^k d^k' 


^-ik-x—ik'-x 


(27r)® + ieko) {k'‘^ + iek^y) 


k+ Du - k^ D++ - ^ ^ ^ 


o I 


and 


(36) 


f d^k d^k' 

^-ik-x—ik'-x 

rip , 1 

~ J (27r)8 

{k‘^ -|- ieko) (A:'2 -|- Ac/cq) 

— 2 — Dii 4“ 2k— k_ -|- ki kj Dij 


X i^5(i + *:') 

Substituting Eqs. dsni) and (jSID into Eq. (ISl yields 


( 37 ) 


d'kd'k' 


^-ik-x—ik'-x 


e = 


(27r)® + ieko) {k'‘^ + iek'^) 


2 U+y 


1 Kk'^--e 


Du + 


(38) 

Since the tensor structure of the correlators D^^ from Eq. m is symmetric under k ^ k', and 
using the last relation in Eq. (Id2j) . without any loss of generality one can write 

D++ = k+k'^f2{k\k'\k-k',kT), (39) 

where f 2 {k\k'\k- k'^kr) is some unknown boost-invariant function, which, due to rapidity 
independence of the problem, depends only on k"^, k'"^, k ■ k' and on the magnitude of the 
transverse momentum kT- In general, dependence of f 2 on k ■ k' might lead to rapidity depen¬ 
dence: however, as we will see below the resulting leading energy density is still boost invariant. 
f2{k\k'\k-k\kT) is symmetric under the interchange /c'^. Similarly 

Du = h{k\k'\k-k\kT) ( 40 ) 


2 A:' - - f — 1 A;2 

- 2 UJ - 


D++ + 


4 U+i U_ k'_l 


and 

kikjDij = f5{k‘^,k'‘^,k-k’,kT) ( 41 ) 

with /i and /s also some symmetric functions under k"^ k'"^. Using these redehnitions we can 

rewrite Eq. dHEl) as 


e 


d^kd^D ^-ik-x-ik'-x 

(27r)® (fc2 -|- ieko) (A)'^ + iek'^) 


1 

2 




h{k\k'\k-k\kT) + 


2 k. 




k+ k'j^ /2(A:^, k'‘^,k- k', kT) + 
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+ 


1 

4 




+ 1 


k',kT) 




5{k + k!)- 


(42) 


For reasons which will become apparent in a moment, we are interested in determining the 
following combination of /’s 


fi{k^ = 0, k'"^ = 0, A; ■ fc' = 0, kr) - = 0, k'"^ = 0, A; ■ A;' = 0, kr) 

2 




/3(A;2 = 0, A;'2 = 0, A; ■ A;'= 0, A;t). 


(43) 


To calcnlate it we compare Eq. (jSED with 


= 


d^kd^k' 


^-ik-x—ik'-x 


[k+rp{k) - kpJi{k)] 


(27r)® (A:^ + ieko) {k'"^ + iek'o) 

X [k\j;(k') - k'^jnk')]y^''^s(k + eh 

which follows from Eq. (Ei. Eqnating the integrands of Eqs. dSHD and we derive 
A:+ k'j^ fi{k'^, k''^, k ■ k\ kx) — k^ k'_^_ f2{k'^, k''^, k ■ k', kx) — 


(44) 


“5 (r + ^) Mk^k'^,k-k',kT) = (45) 

Pntting k = —k' and k‘^ = A;'^ = 0 in Eq. (03) and employing the fact that in covariant gauge 
kPJp{k) = 0 we obtain 


/i(A;^ = 0, k'^ = 0, A: ■ A:' = 0, kx) — kxf 2 {k‘^ = 0, k'"^ = 0, k ■ k' = 0, kx) — 
--sf3(k'‘ = 0,k‘^ = 0,k-k' = 0,kT) = 


(46) 


Finally, since in order to construct the amplitude out of the field given by Eq. (OH|l one needs to 
truncate the held and put the outgoing gluon’s momentum on the mass shell, k"^ = 0, we see that 
J°‘p{k) at A:^ = 0 is nothing but a production amplitude for a real gluon carrying momentum k 
(without convolution with the polarization vector). The corresponding multiplicity distribution 
of the produced gluons is given by 


dN 
d'^k dy 


2(27r)= 


r^{k) r(-k) 


A:2=0 


(47) 


Therefore, 


fi{k‘^ = 0, A:'^ = 0, A: ■ A:' = 0, kx) — kxf 2 {k‘^ = 0, A:'^ = 0, k ■ k' = 0, kx) 


k^ 


fsik^ = 0, k'^ = 0,k-k' = 0, kx) = -2(27r)'’ 


dN 
d‘^k dy 


(48) 
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and, for a cylindrical nnclens, 


1 


= 0, k'^ = 0, ■ fc' = 0, kr) - 4/2(fc = 0, k'^ = 0, fc ■ fc' = 0, kr) 


k^ 


= 0,k'^ = 0,k-k' = 0,kT) 


= -2(27r) 


dN 


d‘^kdyd% 


Now let ns get back to Eq. (P|). Rewriting for each of the /’s 

k ■ k\ kx) = fi{k‘^ = 0, = 0, fc ■ A;' = 0, fcr) + 


(49) 


+ [fi{k‘^, k'‘^, k ■ k\ kx) — fi{k‘^ = 0, = 0, fc ■ A;' = 0, fc^)] (50) 

and keeping only the fi{k‘^ = 0, k'^ = 0, k-k' = 0, kr) in Eq. il we can perform the longitndinal 
momenta (A;+, k-, k'_^_, k'_) integrations with the help of Eq. ()R7jl from Appendix B obtaining 


e ~ J {[Jl{kTT)f + [MkTT)f} 


X 


fi{k^ = 0, k'^ = 0, k ■ k' = 0, kT) — k^f2{k; = 0, k 


/2 


0, k ■ k' = 0, kx) 


~ 7^ = 0, A: ■ A:' = 0, kr) , 

fv'j-' 

which, after employing Eq. (gni) becomes 


(51) 


(Again we have nsed the rapidity-independence of the glnon spectrnm to replace y with 

y.) 

One might worry that the fnnctions fi{k‘^, k''^, k-k', kr) may not have a hnite k'^,k'‘^,k-k' —0 
limit, which would be dangerous for the decomposition of Eq. piji [zg. However, let us remind 
the reader that the quantity J“(A:) dehned in Eq. has the meaning of (truncated) gluon 
production amplitude for the off-shell gluon with virtuality A;^. In the A;^ ^ 0 limit J))(A:) 
becomes the gluon production amplitude for an on-shell gluon, and is indeed hnite. Therefore, 
the correlation functions from Eq. BD, which in the k‘^,k'‘^,k-k' —> 0 limit have the meaning 
of the gluon production amplitude squared (but without the Lorentz index contraction), are 
also hnite in this limit. This implies that the functions /i(A:^, k'"^, k ■ k', kT), dehned in terms of 
various components of in Eqs. dsni), (SOI) and (ED, are hnite in the k‘^,k'‘^,k ■ k' ^ 0 limit. 

For the proper time r much larger than 1/(A:^), with kT the typical transverse momentum 
in the distribution Eq. (IKgi becomes 


e 


T>l/{fcT) 



J ^ d^kdyd% 


1 dE/T 
T dr]d?b' 


(53) 


i.e., it falls oh as 1/r. 
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Therefore, we have shown that the energy density e of a glnon held prodnced in a heavy 
ion collision always has a non-zero term scaling as ~ 1/r. However, to demonstrate that this 
term dominates at late times, we still need to prove that it does not get canceled by the terms 
we left out in writing down the decomposition of Eq. (isnD and keeping the hrst terms only. 
Thus we have to analyze the contribution arising from substituting the terms from the square 
brackets of Eq. dsni) into Eq. (I42|l . We need to show that such contributions fall off faster than 
1/r, and, therefore, can be neglected at late times. Here we will demonstrate that this is true 
for one of the terms — the /i-term in Eq. (I42j) . The proof for the other two terms on the right 
hand side of Eq. (P|) would be analogous. 

Substituting the square brackets from Eq. dSHD into Eq. dll we obtain the following con¬ 
tribution to the energy density, which we have to prove to be small; 

1 r dk'j_ dk'_ ^-ik-x-ik'-x 

2 J (27r)® S± {k? -f ieko) -|- iek'o) 



X [/i(A;^, k ■ k', kx) — fi{k‘^ = 0, k'"^ = 0, k ■ k' = 0, kx)]- (54) 

For a wide range of amplitudes one can write 

/i(/c^, k'‘^, k ■ k', kx) — fi{k‘^ = 0, k'"^ = 0,k ■ k' = 0, kx) = {k"^ k''^)^^ k'‘^, k ■ k\ kx) + 

+ [{k + k'ff^ g^^\k\ k'\ k ■ k\ kx), (55) 

where g'd'){k‘^ = 0 , = 0 , A; ■ fc' = 0 , /ct) ^ 0 , g'^‘^\k‘^ = 0 , = 0 , /c ■ fc' = 0 , /ct) ^ 0 , and 

Ai,A2 > 0 . In arriving at Eq. dsi we have also used the fact that fi{k‘^,k'‘^,k ■ k',kx) = 
/i(fc'^, k ■ k', kx), which follows from the k k' symmetry in the dehnition of /i(A;^, k''^, k ■ 
k', kx) given by Eq. dm along with Eq. m- In Eq. dm we put a power of {k + k')^ instead 
of a power of fc ■ A;' in front of g^‘^\ Similarly to Eq. dm we write 

g^^\k\ k'^, k ■ k', kx) = g^"\k^ = 0, k'^ = 0,k-k' = 0, kx)+ 

+ [g^^\k‘^, A:'^, k ■ k', kx) — g^^\k'^ = 0, k'"^ = 0, A; ■ A:' = 0, kx)] (56) 

for i = 1,2. Substituting the first term on the right hand side of Eq. dm for g^^'> into the 
hrst term on the right hand side of Eq. dm, and then plugging the resulting contribution into 
Eq. (im yields 


1 

2 


d‘^k dk'j^ dk'_ e 

(27r)6 S'j_ (A;2 -|- iek^) (A;'^ -|- iek'^) 



(/•'•(O.O.O.fc'r). 

(57) 


Performing the A:_|_, A:_, A)(,_, k'_ integrations in Eq. dm with the help of Eq. (EH) in Appendix A 
and Eq. (IB5jl with A = 1 in Appendix B we obtain 


g27riAi 

8^xr(i-Ai)2 


d?k 


g^^\0,0,0,kx) kx 



which, for r S> 1/ (kx), scales as 


1 

7T+2a7’ 


(58) 


(59) 
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and is, therefore, negligibly small at late proper times compared to the leading contribntion 
to energy density given by Eq. (Here we assnme that the typical transverse momentnm 

(fcr) is the same for in Eq. (j5^ and for g^^\0,0,kT) in Eq. (j58|l ; both fnnctions 

resnlt from expanding the same amplitnde in powers of k''^, and no new scale can arise from 
snch an expansion, which jnstihes onr assnmption.) The particnlar way of regnlarizing the k‘^ 
branch cnt nsed in Eq. m and Eq. ED is not essential for arriving at Eq. dSHD, since other 
regnlarizations wonld yield the same resnlt. 

The second term on the right hand side of Eq. (jSi gives a similarly small contribntion. To 
see this we snbstitnte the hrst term on the right hand side of Eq. (jsni) for into the second 
term on the right hand side of Eq. dSHD, and then snbstitnte the resnlt into Eq. (jSl obtaining 


d^kdkk dk' 


-.—ik-x—ik'-x 


(27r)®2S'_L {k"^ + ieko) {k'"^ + iek'o] 






1 r d^k 

= J {[Ji{kTr)f +[MkTT)f} ■ (60) 

For r S> l/(/cr) the integral on the right of Eq. (IHn|l scales as ~ 1/r: applying the derivatives 
we see that the whole expression in Eq. dsni) scales as 


1 

T-1+2A2 ’ 


( 61 ) 


and is also negligibly small at late proper times compared to the leading contribntion to energy 
density given by Eq. 

For the second term on the right hand side of Eq. (j56|l . gk\k^, /c'^, k-k', kx) —g^^\0, 0,0, /cr) 
with i = l(or i = 2), one can repeat the procednre ontlined above for /i(fc^, k'^, k ■ k', kr) — 
/i(0, 0, 0, /ct), nsing the redehnition jnst like in Eq. and showing that the leading term in the 
resnlting decomposition, similar to Eq. dsni), falls off faster with r than Eq. (Ei (or Eq. (ED))- 
Iterating the procednre wonld generate a series of corrections falling off at progressively higher 
powers of 1/r, all of which could be neglected at r 1/(fcr)- 

Of course the assumption of Eq. dSHD, while quite general, does not include all the possibil¬ 
ities. One might imagine other ways for fi{k^, A;'^, k ■ k', kr) — /i(0, 0, 0, kr) to approach zero 
as fc ■ fc' —> 0: it might scale as l/(lnA:^ In or, less likely, as 

case, Eq. (ED suggests that each power of k"^ (or each power of k'"^ or of A: ■ k') gives a power 
of 1/r for energy density e in coordinate space: k"^ —> 1/r. (Indeed the powers of kT do not 
modify the r-dependence of at all.) Therefore, one may argue that after the momentum 
integration is done in Eq. (IS2D, the /i(A;^, A:'^, k ■ A:', kx) — /i(0, 0, 0, kx) term, when substituted 
into Eq. (I42|l . yields approximately the following contribution to e 


1 

r 


fi 


111 


T T T 


5 (kr) 


/i(0,0,0, {kr)) 


(62) 


which falls off faster than 1/r and can thus be neglected compared to Eq. (jbdj) . This conclusion 
is natural, since the term in Eq. or, equivalently, the second term on the right hand side 

of Eq. (IHn|l . does not contribute to the production cross section, as follows from Eq. which 
is another way of saying that it is not important at late times. 
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The proofs that the contributions to energy density e generated by substituting / 2 (fc^, k'"^, k- 
fc', kx) —/2(0 , 0, 0, /ct) and fsik'^, k''^, k ■ k', kx) —/3(0, 0, 0, kx) instead of /2 and /s into Eq. (11^ 
are also subleading at large r can be constructed by analogy to the above. 

Finally, we have to comment on our use of the Abelian part of only in Eq. (I2ni) and 
throughout this Section. Including the non-Abelian parts of the held strength tensor would 
generate higher powers of A“ in the dehnition of Using Eq. (1^ those extra powers 
can be rewritten as extra integrals over k" and k'" in the extra terms which would be added 
to Eq. (1^ . Due to Eq. (EH), each of these extra integrals would (at least) generate a Bessel 
function J^/^^kxT), which at large r scales as (I/a/t) cos^krT -|- |A — ^). Even without the 
cosine, one can immediately see that the cubic in A“ term in would fall off at least like 
at large r. The quadric terms would fall off at least like l/n^. Both of these terms would 
be negligibly small compared to the leading quadratic term scaling as 1/r shown in Eq. (jh.'fjl . 

Eq. has a straightforward physical interpretation. Every Feynman diagram has a hnal 
state in which the particles are propagating as non-interacting plane waves until the inhnite 
late times. Indeed the energy density of such a “free-streaming” state scales as ~ 1/r, and this 
is exactly what Eq. represents. 

Therefore, in this Section we have proven that in the rapidity-independent case, dehned by 
Eqs. m and m for the total rapidity interval in the collision of two very large nuclei, at any 
order in the perturbative expansion in the strong coupling g, the resulting gluon held’s energy 
density has a non-vanishing term which is dominant at late times giving e ~ 1/r flhdj) . Hence 
it appears that, in this boost-invariant case, thermalization leading to Bjorken hydrodynamic 
description of the evolution of produced gluonic system, can not result from resummation of 
perturbative QCD diagrams. 

4 Generalization to the Rapidity-Dependent Case 

For rapidity intervals U > — in heavy ion collisions the quantum evolution corrections IZniEBl 
EH ESI would become important making the produced particle distribution rapidity dependent. 
Below we are first going to argue that rapidity-dependent hydrodynamic description may only 
change the e ~ l/r^^^ scaling of the ideal Bjorken energy density to a higher power, e ~ 
l/.^4/3+A_ yy-g then demonstrate that the rapidity-dependent quantum corrections, such as 
the ones introduced by the BFKL evolution eg, would not modify the e ~ 1/r scaling derived 
in the previous Section. 

4.1 Rapidity-Dependent Hydrodynamics 

In the rapidity dependent case the most general form of the energy-momentum tensor is given 
by the equation similar to Eq. (USD 

= A{t, 7 ]) + B{t, 7 ]) + u„v^) + C{t, 7 ]) -f D{t, 7 ]) (63) 

with and still given by Eqs. dn and dUD and where now all the coefficients A, B,C, D 
are also functions of the space-time rapidity t]. Due to this 77 -dependence the -|- — symmetry 

argument no longer applies in general. However, it still holds at mid-rapidity (77 = 0) for a 
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collision of two identical nuclei leading to 


B(r, ?7 = 0) = 0. (64) 

Applying the conservation of energy-momentum tensor condition (|22D to the tensor in Eq. (IMD 
yields 


TdrB - 2dr,D + 2d^C + 2B = 0 
2TdrA + dr^B + 2TdrD + 2A + 2C = 0. (65) 

The energy-momentum tensor in Eq. would describe a hydrodynamic system if it could be 
reduced to the standard hydrodynamic form 

= {e + p)w^Wy-pg^u, ( 66 ) 

where is the four-vector of the fluid velocity, = 1. For the 1-|-1-dimensional expansion 

of the system created in a collisions of two very large nuclei considered here the fluid velocity 
has zero transverse component, w = 0, such that = (w+,tc_,0). Matching Eq. onto 
Eq. (|F)!)|l we obtain 

A = e + p + C (67) 

and 

D = -p. (68) 

For the hydrodynamic energy momentum tensor (IHH]) the following relation holds 

T++T__ = (T+_+p)2, (69) 


leading to a constraint 



(70) 


Combining Eq. (inni) and Eq. (HOD with the equation of state relating e and p, would give us 
a complete set of rapidity-dependent hydrodynamic equations. However the resulting system 
of equations is nonlinear and is hard to solve analytically. Instead we are going to construct 
a perturbative solution for small rapidity-dependent corrections to Bjorken hydrodynamics 
PP. We begin by noting that, since i? = 0 in the boost-invariant case considered in the 
previous Section, we can assume that non-zero B reflects the deviation from the ideal Bjorken 
hydrodynamics, and could be assumed small if we are interested in small corrections to the 
latter. Assuming that B A and keeping only linear in B corrections allows us to neglect C, 
since, due to Eq. dZOD, C ^ B^. Than, using Eqs. (EH) and (EHl) in Eq. (jHil) yields 


r drB + 2 drjP + 2B = 0 
2{Tdre + e + p) + dr,B = 0 . 


(71) 


We are interested in the solution for the ideal gas equation of state: therefore we put e = 3p. 
The most general solution of Eq. m satisfying the condition of Eq. (El and mapping back 
onto Bjorken hydrodynamic behavior for small B is 


e = Co cos 


(VAp) 


(5-v^n^) 


(72) 
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with 


B 


2 VI - 3A - 1 
“3^° Vl 


sin [VAt]) — 

T3 


1 

(5-V1-3A) ’ 


(73) 


where A and eo are arbitrary constants. The corresponding flow velocity components are given 
by 


x± 

w± ~ — 

T 


1 ± 


3B 

Ye 


(74) 


Looking at the solution given by Eq. (|7^ one may wonder why the energy density is not positive 
deflnite. Indeed for A < 0 the energy density e from Eq. (izg) becomes positive deflnite, since 
cos{y/Ari) would be replaced by cosh(y|A| rj). However, the resulting rapidity distribution of 
energy density would increase as one moves further away from mid-rapidity, which is unphysical. 
Therefore one has to have A > 0. Resolution of the positivity problem for e comes from the 
necessity to satisfy the i? -C A assumption which we have made at the beginning of this 
calculation. It translates into i? <C e condition, which is satisfied by Eqs. and dZSD only 
if Va// 1. Since in this Section we are interested in large rapidity intervals, rj y>l/a., 
the a/A ?7 -C 1 requires that A < a* -C 1. Hence, for large rapidities, Eqs. (ff^ and are 
valid only at the lowest order in A 


e 




(75) 


B ^ (76) 

where we did not expand since, at late times, A Inr does not have to be small for R -C e 

condition to hold. For small VA r] the energy density in Eq. ea is indeed positive. 

Eq. ()75fl has an important feature which we would like to emphasize: since A > 0, it 
shows that the energy density of the boost-non-invariant ideal hydrodynamics falls off with r 
faster than the energy density of the boost-invariant ideal Bjorken hydrodynamics pp. Here 
we have proven it only for a small rapidity-dependent perturbation of the Bjorken solution. 
However one should expect our conclusion to hold in a general case of a rapidity-dependent 
hydrodynamics. In the case of a rapidity-dependent hydrodynamics, the longitudinal pressure 
is higher than in the boost-invariant Bjorken case, generating the longitudinal acceleration of 
the flow (see e.g. Eq. (d)- The central-rapidity high-density system starts expanding faster 
than in Bjorken case, leading to a faster depletion of the energy density with r. In other words, 
once the longitudinal homogeneity of pure Bjorken hydrodynamics is broken by some rapidity- 
dependent phenomena, the system starts doing more work in the longitudinal direction than 
it was doing in pure Bjorken hydrodynamics case, and this leads to a faster decrease of energy 
density with proper time.^ 


4.2 Rapidity-Dependent Energy Density 

Here we are going to generalize the argument of Sect. 13.31 to the case of rapidity-dependent 
gluon fields. It is impossible to define co-moving energy density for a general energy-momentum 
tensor like the one given in Eq. dnsi), since, in the general not necessarily hydrodynamic case, 

^The author would like to thank Ulrich Heinz for explaining to him this argument. 
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one can not define the co-moving frame, and in the case of hydrodynamics (IM|) one needs to 
know the flow velocity to define the co-moving frame, which is impossible to do withont solving 
the hydrodynamics equations Therefore we will restrict our analysis to the case of mid¬ 

rapidity, rj = 0, where, for a collision of two identical nuclei, the co-moving frame is just the 
center of mass frame of the two nuclei. There Eq. (jSSI) would apply, such that 


e{T,T] = 0) = ^ T++(r,r 7 = 0)+T+_(r,r/= 0) 

Repeating the steps from Section l!01 we write 


T++{T,r] = 0)+T+_{T,r] = 0). (77) 


e(T,rj = 0) = y 


d^kd^k' 


^-ik-x—ik'-x 


(27r)® (k‘^ + ieko) (k'"^ + iek'^) 


r]=0 


1 

2 




X fi(k‘^, k+, k'‘^, k'_^_, k ■ k', kx) -|- 


2k_ k'_ 


1 

2 



k+ k'^ f2(k‘^, k+, k''^, k'_^_, k ■ k', kx) 


+ 


1 

4 




h(k\k+,k'\k'^,k-k\kx) 



d(k + ^). 


(78) 


where now, in the rapidity dependent case, /ds are functions of k± and k'^ as well. However, 
since we can always rewrite /c_ = -|- k‘^)/2k+ and k'_ = (k'^ + k‘^)/2k'^, we put only k^ and 

k'_^_ in the arguments of the functions /*. 

Rapidity-dependent quantum evolution corrections come in as logarithms of Bjorken x vari¬ 
able If p+ is a large longitudinal momentum carried by a nucleon in the nucleus moving 
in the -I— direction, than x = The rapidity-dependent corrections would then bring in 

powers of cts In l/x = lnp+/fc+. Resummation of such corrections for the gluon production 
amplitudes generates powers of 1/a:, or, equivalently, p+/k+. Therefore, to verify whether such 
corrections modify the r-dependence of e, we can consider the following general form for the 
functions //s 


fi{k‘^,k+,k'‘^,k'^,k ■ k',kx) 




fi(k‘^,k'‘^,k- k',kx), 


(79) 


where we again used the fact that /’s are symmetric under k k' interchange. For simplicity 
we assume the power A to be the same for /i, /2 and /s: this assumption is not crucial and can 
be easily relaxed. The logarithmic corrections to //s, i.e., terms with lnp_|_/fc_|_ and lnp_|_/fc(_, 
can be obtained from //s in Eq. (17^ by differentiating it with respect to A. Indeed that would 
give only logarithms like ln(p^/fc+A;(_), but not lrLk^/k'_^_: while we are quite confident that the 
latter terms never appear in perturbation theory, our approach can be easily generalized to 
include both types of logarithms by putting different powers for / k^ and / k'_^_ factors in 
Eq. ()79j) . (A careful reader may worry that the p+-dependence was never explicitly shown in 
the argument of //s and was explicitly assumed there: in fact, due to boost-invariance, the 
p+-dependence enters in //s only through the ratios of p+/fc+ and p+/k'_^_. In the rapidity- 
independent case of Section l!01 //s were independent of p+, which corresponds to the eikonal 
limit.) 
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Substituting /j’s from Eq. dZi into Eq. dZHD we obtain 


e{T,r] = 0) = y 


d'^kd^k' 


^-ik-x—ik'-x 


( 27 r)® {k? + ieko) {k'‘^ + iek'^) 


r)=0 


1 

2 




X fi{k'^, k ■ k', kx) + 


2 k_ k'_ 


1 

2 



k+k'^ f2{k‘^,k'‘^,k ■ k',kT) + 


+ 


1 

4 




f3{k^,k'^,k- k\kT) 


U+ fcVyl 



5(fc + fc')- 


(80) 


To perform the longitudinal momentum integrals we will use the integral in Eqs. (IBljl and (IB6f) 
of Appendix B. There we can see that, similar to the rapidity-independent case, each positive 
extra power of (or k'^ or k ■ k') gives a power of 1/r. Therefore, we again are interested in 
contribution of fi{k^ = 0, A:+, = 0, k'j^, k ■ k' = 0, kx) as in the terms giving the leading-r 

behavior. Similar to Eq. ()40j) we can write 


(P+ P+ 


fc I — fci 


/l( 0 , 0 , 0 , kx) - 4 / 2 ( 0 , 0 , 0 , kx) - -rj / 3 ( 0 , 0 , 0 , kx) 

rCj-’ 


= -2{2Tif 


dN 


d‘^k dy d% 
(81) 

which shows that this combination of //s is not zero. Using Eq. m in Eq. dHOl), where we put 
k^ = = 0 in the arguments of all //s, and performing the longitudinal integrations using 

the formulas from Appendix B yields for the leading term in energy density 


e(r, 77 = 0) 


71 

2 



dN 

d?kdr]d% 




(82) 


In arriving at Eq. (IH^ we have noticed that, according to Eq. (jBOj, each power of /c+ or fc/ 
gives a power of kx eP j V2 after the integration. For on-mass shell gluons in Eq. m one has 
k+ = kxe^/^/2. Therefore, the powers of kxe^/\/2 were absorbed in ^y just replacing 

y ^T]. 

The large-r asymptotics of Eq. (IH21) is the same as in Eq. (jS3I): 


e(r,r/ = 0) 


T>l/(fcT> 


T 



dN 

d‘^k dr] d% 



1 dEx 

T dyd^b 


r]=0 


(83) 


Therefore, we have proven that even in the rapidity-dependent case the mid-rapidity energy 
density given by the Feynman diagrams falls off as 1/r at large r. This conclusion could be 
easily derived by just analyzing Eqs. EO and dHEl): one can see there that each extra power 
of A:+ does not bring in any new powers of r, and only modifies the order of the Bessel function, 
which can not change the r-dependence, as follows from Eq. (IH^ . 

Eq. shows that the scaling of energy density of the rapidity-dependent solution of 
the hydrodynamics equations given by Eq. dZSl), e ~ l/r3+2, can not come from the leading 
contribution of Feynman diagrams. Indeed, the subleading contributions may still lead to 
energy density falling off with r faster than 1/r: however, due to Eq. (lAlll . such contributions 
must come in with extra positive powers of k'^ in momentum space. They would go to zero in 
the on-mass shell —> 0 limit and, thus, would not contribute to the production cross section. 
Such corrections are probably irrelevant for all physical observables. 
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5 Including Quarks 


To generalize our conclusion to massless quarks we will restrict our discussion to rapidity- 
independent case only: generalization to the rapidity-dependent case can be easily done fol¬ 
lowing the procedure outlined in Section 14.21 We start with the energy-momentum tensor for 
a single massless quark flavor: 




(84) 


The corresponding energy density is given by Eq. O, which we again want to rewrite as a 
double integral, just like Eq. ()42ll . by Fourier transforming the quark held 


and 


%Ij{x) = i 


= i 


' d'^k /c ■ 7 

(27r)4 P + teko "" ^ 

(85) 

d'^k' r,,,. A;'■ 7 

( 86 ) 


with ^{k) and ^{k') some spinors. In the following, similar to Section ESI we will keep only the 
Abelian part of The non-Abelian corrections are suppressed at late times and can be 

neglected, since, just like in Section 1,01 they fall off faster than the Abelian term by at least a 
factor of 1 /a/t. Replacing the covariant derivatives in Eq. by a regular derivative 
and substituting Eqs. and (jHED in it we obtain 


rjiquark _ __ 

- 2 


d^kd'^k' 


^-ik-x—ik'-x 


l{k') k' ■ 7 ( 7 ^ ky + 7 ^ kfj) k ■ 1 ^{k)) . (87) 


(27r)® -I- ze/co) {k'"^ -|- ze/cg) 

Rewriting 

lik') k' ■ 77/4^ ■ 7C(^))) = - - k'^) k ■ k', kr) + + k'^) h2{k‘^, k ■ k\ kr) 

in Eq. (|H7|) and using Eq. we write 

d^kd^k' 


^-ik-x—ik'-x 


^quarkf^\ _ 


X 


\ (^) ~ k+ + k+k- - k- + k+ k'_) 


{2'kY {kP' -|- ieko) [k'"^ -|- zefcg) 

h{kYk'Yk-k\kT) + 


\ (“) {k++k'_^) k++k+k_ + ^ {k'_^_ k-+k+k'_) 


h2{kYk'‘^,k-k',kT) 


S± 


6{k+k!). (89) 


Similar to Section EiSl by putting k = —k' and k"^ = = 0 in Eq. 

^ hYk^ = 0, k'^ = 0, fc ■ A:' = 0, fcr) = ^ Ui{-k) k -^ak) 

o \ O I ' ^ 


we derive 


dm 


A;2=o ^ ^ d‘^kdyd%^ 


(90) 
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where ^2kdy(Pb multiplicity of the produced quarks. Arguing, just like we did for gluons, 

that the leading-r behavior for the quark energy density is given by hi(/c^ = 0, k'"^ = 0, k ■ k' = 
0, kr) and h 2 {k‘^ = 0, = 0, k ■ k' = 0, kx) in Eq. (j89|l we can integrate over the longitudinal 

momenta in Eq. (iHnD using Eq. dEzj) in Appendix B obtaining 

^ ^ j (fk {-Mf^TT) J2{kTT) + 2[Ji(/Ctt)]^ + [Jo{kTr)f} + 

+ J ^ 2 ( 0 , 0, 0, kr) kl J2{kTT) - 2[Ji{kTT)]^ + [MkrT)]^} , (91) 

where we also substituted space-time rapidity rj instead of y in the quark multiplicity distribu¬ 
tion, which makes no difference in the boost invariant case we consider. We can assume that 
h 2 {k‘^ = 0, k'"^ = 0, k ■ k' = 0, kr) is a steeply falling function of kx for kx 3> {kx) ~ Qs- The 
assumption is justihed since h 2 {k‘^ = 0, = 0, fc ■ A;' = 0, kx) comes from the same amplitude 

that gave hi(A;^ = 0, k'"^ = 0, /c ■ fc' = 0, kx), which is equal to the quark spectrum, as shown in 
Eq. (inn|l , which in turn is always a steeply falling function of kx scaling at least like ~ 1 / for 
kx above some scale {kx) ~ Qs- By the same argument h 2 {k‘^ = 0, k''^ = 0, /c ■ /c' = 0, kx) should 
be regular (or at most logarithmically divergent) at kx = 0. Using these assumptions we can 
argue that the integral in the second term on the right hand side of Eq. (EH) is dominated by 
kx ^ Qs, which allows us to rewrite it as 

—^2(0,0,0, Q^) [ dkxkx {-JoikxT) J2{kxT)-2[Ji{kxT)f+ [MkxT)f} ■ (92) 

lO TT JO '' 

Performing the integration in Eq. dni one would obtain a linear combination of hypergeometric 
functions, which can be shown to fall off at least as ~ I/'t^ at large r. Therefore the second 
term on the right hand side of Eq. dnH) falls off with r at least as ~ 1 /t^, and can be neglected 
if the hrst term on the right hand side of Eq. (EH) falls off with r slower than l/rT 

This can be easily verihed. At large r the hrst term on the right hand side of Eq. EH) gives 


quark, 


^) 


- fd^k 


dm 


r>l/(fcT> 


r 


d?kdr]d% 


kx — — 


1 


T dr] d% ’ 


(93) 


since the factor of 2 accounts for the anti-quark contribution. We can see that indeed the term 
in Eq. El is negligibly small compared to Eq. El, which gives us the dominant contribution 
to the energy density due to quarks at late times r. (Of course the typical transverse momentum 
{kx) in Eq. (El does not have to be exactly equal to the similar typical momentum for gluons 
in Eq. El: however, the difference between the two is usually given by the ratio of the Casimir 
operators, which is just a constant (4/9) and does not change our argument above.) 

We have shown that inclusion of massless quarks does not change the conclusion of the 
previous Sections that the leading diagrammatic contribution to energy density scales as 1/r 
at large r for any order in the coupling g. Therefore, it appears that inclusion of quarks does 
not affect the onset of thermalization. 


6 Conclusions 

We have shown above in Eqs. (El and El that gluon and quark helds generated by Feynman 
diagrams in high energy heavy ion collisions lead to energy densities scaling as e ~ 1/r at 
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r ^ l/{kj)- In the saturation/Color Glass picture of heavy ion collisions, the typical mo¬ 
mentum (/c't) is proportional to the saturation scale Qs- Therefore, the 1/r scaling of energy 
density sets in at relatively early times r ~ 1/Qs (even though throughout the paper we called 
these proper times late times). The remaining evolution of the system in the perturbative 
scenario considered above, characterized by e ~ 1/r scaling, is reminiscent of the so-called free 
streaming, where the system simply falls apart without particles interacting. 

To explain how this happens, let us provide a diagrammatic interpretation of our conclusion 
of e ~ 1/r scaling. Let us imagine a general gluon held produced in a heavy ion collision as 
shown in Fig.|2K. There the gluon held is hrst produced in the nuclear collision at r = 0 denoted 
by the ® sign. The cross at the other end denotes the later point in r where we measure the 
energy density of the gluon held. In the evolution of the system the gluon interacts with other 
gluon helds produced in the collision in all possible ways, as shown in Fig. |2K. However, the 
proper times of these interactions are not hxed: they are integrated over the whole range of 
r. Interactions may also happen at diherent impact parameters, which are also integrated 
over. The e ~ 1/r scaling conclusion from Eqs. and appears to indicate that the 

dominant diagrams are given by Fig.EB, where all the interactions happen at early times, after 
which the system simply falls apart. In other words, the integrations over proper times of the 
interactions in Fig. |21^ are dominated by early times of Fig. |2l3. (k-r), or Qg, being the only 
scale in the problem, sets the typical time scale for the end of interaction period and the onset 
of free streaming, tq ~ 1/Qs- Such behavior has been previously observed in the numerical 
simulations of the classical gluon helds isziEni. 

Another way to physically understand our conclusion of e ~ 1/r scaling is as follows. At any 
order in the coupling constant the gluon (or quark) held has a diagram (or several diagrams) 
which is (are) non-zero if the gluon is put on mass-shell. This is just a statement that gluon 
(or quark) multiplicity distribution can be expanded in a perturbation series in As we have 
shown above in Sections lOl R3I and El such diagrams always give energy density scaling as 
1/r. Each diagram is dominated by the on-shell particles free streaming away, which always 
leads to e ~ 1/r. 

Therefore we have shown that the onset of thermalization and the subsequent Bjorken or 
rapidity-dependent hydrodynamic expansion of the system of quarks and gluons produced in 




A B 

Figure 2: (A) Gluon held produced in a collision with all the interactions throughout its 
proper time evolution. (B) The dominant contribution to the gluon held comes from early-time 
interactions. 
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heavy ion collisions can not result from summation of Feynman diagrams. Nevertheless, there 
exists a solid phenomenological evidence for the strong hnal state interactions 01110131101 
ITU IT 2 I ITdl HH im Cni im CH] and hydrodynamic behavior |2l 0 IH 0 of the system produced 
in heavy ion collisions at RHIC, indicating a formation of strongly interacting quark-gluon 
plasma (QGP). To reconcile it with the above argument that Feynman diagrams lead to a free 
streaming behavior for both quarks and gluons, one must conclude that non-perturbative QCD 
effects are instrumental in QGP formation at RHIG. These could be the the non-perturbative 
effects associated with infrared modes having momenta of the order of Aqqd which can not 
be represented by Feynman diagrams. Therefore, the above argument does not apply to such 
modes. Alternatively, the non-perturbative effects might be of the nature similar to the ultra- 
soft modes in hnite temperature non-Abelian field theories, which have momenta of the order 
of T with T the temperature of the system. It is well-known that resummation of ultra-soft 
modes is a non-perturbative problem in finite temperature QGD eg. It is also known that 
ultra-soft modes are very important for many physical observables for equilibrium QGD matter 
at hnite temperature [zniiHiiini. If they are important for equilibrium QGD matter, it would 
be natural to suggest that the ultra-soft modes could also play a major role in non-equilibrium 
phenomena such as the onset of thermalization. However, a more careful analysis of the issue 
is needed in order to draw any conclusions. Such analysis is beyond the scope of this paper. 

Distinguishing which one of the two types of non-perturbative effects plays a more important 
role in the process of thermalization would also be important for our understanding of LHG 
heavy ion data. The non-perturbative effects characterized by the scale Aqcd are likely to 
be of little importance at LHG where the saturation scale Qs is predicted to be much larger 
than Aqcd shifting most partons away from the infrared region. At the same time, the non- 
perturbative ultra-soft modes carrying momenta T may remain important even at high LHG 
energies if the relevant temperature scales with the saturation scale, T ~ increasing at high 
energy. 
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Appendix A 

Here we are going to prove the following formula 

/ OO 

dk+ dk_ + ieko)^-^ -- 

-OO 


27r2 


'2k7 


F(l-A) V T 


e'-^J_A(fcTr) (Al) 
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with = 2 ^ and for x+ > 0, x_ > 0, and A > 0. Let us first rewrite the integral 

flAlfl as 


A_i dk+dk_e / 

J-oc {k+ + iey~^ y 

Dehning k- = k_ — /2k^ we write 

dk^ 


k^ 


A-l 


2(/c+ + ie) 


-T + 


1 = 2 


A-l 


oo (A^-|- “t“ 


-A 


e * 24 +-^+ r dk-e-^'^-^+{k-+ie)^-\ 


(A2) 


(A3) 


The k- integral can be easily performed by distorting the integration contour around the branch 
cut. We obtain 


_2A-. ^ 


d/c I 




(A4) 


r(l-A) + 7-00 (fc+ + ie)i-^ 

Expanding the second term in the power of the exponent in Eq. (IA4I1 in a Taylor series we write 




“ 1 f-ik^x+Y dk^ 


r(l-A) 


^—ikj^X- 


n=0 


n\ 


1—00 {k-\- T fe) 


— A 


(A5) 


Performing the /c+ integration just like we did the fc_ integral above yields 

j ^ -pA-i _1_ / -k x+x_ 

r(l-A)^+-^ io^!r(n + i-A) 1, 2 , 


(A6) 


Remembering that 2x+x- = and performing the summation over n we obtain Eq. m as 
desired. 


Appendix B 

Our goal in this appendix is to perform the following integration 

/ CXD 

dk+ dk_ e-*^+^— {k^ + iekoY~^ {k+ + ie)\ 

-OO 

Repeating the steps from Appendix A which led to Eq. (jA4jl we write 

dkx- 


^TTl 2 ^ /*00 

J = -2A-i 5^e2 A f 


—ik+x— —z o ~ , . xjt- 

^ 2/c I +l€ ' 


r(l-A) + 7-00 (A;+ + ie)i-A-A 

Expanding the second term in the exponent yields 


J = -2‘ 


27rie*2 


j-^A 


r(i-A) 
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-A 
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-iY x+ 


dko 


^—ik+x- 


n\ 


' —OO (/c-(- T 


A—A 


Performing the /c_|_-integration we obtain 
(27r)2 ^ 


J = -2" 


P(l-A) 


(x+a;_) ‘^xY E 


—Y x+ X- 


„3i "'r(n + 1 - A - A) 


(Bl) 
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(B3) 


(B4) 
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which, after summing over n gives 
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27r2 

r(i-A) 



kx T 

2x_ 


A 


giTT A+i^A 


J-A-\{kTT). 


Noting that x± = re^^/ a/ 2 with t] the space-time rapidity we rewrite Eq. ()B5jl as 


J 


27r2 

r(i-A) 




gA?7 giTT A+if A 




(B5) 


(B6) 


As one can see from Eq. (|Bi, extra powers of in Eq. CD as opposed to Eq. (EH) do not 
bring in any extra inverse powers of r: they only modify the order of the Bessel function. 

Finally, let us list here another useful integral, which can be easily obtained by direct 
integration 


dk_^_ dk_ ——- - — 

-oo k^ + leko 


klk^ = -2 71^ 


' ikTT^ 
2x_ 


—i kr t' 
, 2x+ ^ 




(B7) 


where n and m are integers. 
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